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Abstract. We give a new algorithm that quickly finds smooth neighbors. 

I. Introduction 

We say that a number is z-smooth if none of its prime factors exceed z. In this paper we 
search for solutions 6 of 

p | 6(6 + 1) p < z. (1) 

If 6 is a solution of (OQ), then we refer to the pair (6,6 + 1) as ^-smooth neighbors. It has 
been known since the work of Stormer in 1898 [S] that for any z there are only finitely many 
^-smooth neighbors. In 1964, D. H. Lehmer [LI] found all 869 of the 41-smooth neighbors. 
To do this he improved upon Stormer's method, which relied on solving a finite number of 
Pell's equations. In fact Lehmer showed that if b(b + 1) is ^-smooth then 46(6 + 1) is the u y 
part" of the nth power of the fundamental solution xq + y/dyo of the Pell's equation 

x 2 - dy 2 = 1 

where d is squarefree and z-smooth, and where n < (d + l)/2. Lehmer solved all of these 
Pell's equations to see which led to 41-smooth neighbors. 

In 2011 Luca and Najman [LN] used a modified version Lehmer's method to find 100- 
smooth neighbors. In a calculation that took 15 days on a quad-core 2.66 GHz processor, they 
found 13,325 neighbors and claimed that this was all of the possible 100-smooth neighbors. 
The calculation was especially difficult because solutions to Pell's equations for squarefree 
100-smooth integers can have as many as 10 10 digits. Consequently they had to use a special 
method to even represent the solutions. In an erratum [LN1], they found 49 more solutions 
that they had missed previously. 

We have found a fast, amazingly simple algorithm that finds almost all z-smooth neighbors 
much more quickly. In fact, when we ran our method to find 100-smooth neighbors, it 
completed in 20 minutes (on a similar machine as Luca and Najman used) and found 13,333 
100-smooth neighbors. We were missing 37 solutions that Luca and Najman found. 

Subsequently we searched for all 200-smooth neighbors. This computation took about 
2 weeks and produced a list of 346,192 solutions. This list included all but one of the 
solutions from the (completed) Luca-Najman list. We determined that this 100-smooth 
number missing from our list would be found using our method when searching for 227- 
smooth neighbors. 
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A plot of the logarithms of the first members from our list of 200-smooth neighbors suggests 
that they are normally distributed. 

2. The algorithm 

Suppose we have a set S of positive integers. For any two elements b < B of S form the 
ratio 

A h B + 1 
J' ~ 6+1 X B 

where gcd(/3, (3') = 1. Sometimes it will be the case that 0' — (3 + 1, for example 

15 _ 3 5 

16 ~ 4 X 4' 

We are particularly interested when this happens, i.e. we are interested in the solutions of 

P + 1 b+ 1 B K ' 

where b and B are in 5*. Given S, we form a new set, S', which is the union of 5* and all of 
the solutions (3 to §Z§. We can repeat this process to form (S')' = S" and so on. Ultimately, 
by Stormer's theorem[S], we will arrive at a set (meaning n iterations of priming) for 
which (S^Y = S( n \ i.e. there are no new solutions to be added. We let 

5(S) 

denote this set that can no longer be enlarged by this process. As an example, suppose 
that 

S = {1,2,3,4,5}. 

Then, it is easy to check that 

S' = {1,2,3,4,5,8,9,15,24}. 

and 

S" = {1,2,3,4,5,8,9,15,24,80}. 
After that we have 5"" = S", so 

£({1,2,3,4,5}) = {1,2,3,4,5,8,9,15,24,80}. 

Recall that Lehmer [LI] gave a complete list of 41-smooth solutions to ([T|). In particular, 
we can see from his list that the above set {1,2,3,4,5,8,9,15,24,80} is the complete list 
of 5-smooth solutions to (CE]). In other words, we found all of the z = 5 solutions to §B) by 
starting with the set {1, 2, 3, 4, 5} and then repeatedly adding in solutions to (T5]). This good 
fortune is not always the case. For example, 

5({1, 2, 3, 4, 5, 6, 7}) = {1,2, 3, 4, 5, 6, 7, 8, 9, 14, 15, 20, 24, 27, 35, 48, 49, 63, 80, 125, 224, 2400} 

whereas from Lehmer's table we see that the complete set of 7-smooth solutions to (OQ) 
includes all of these numbers along with 4374. However, it is the case that 

4374 e {1,2,3,4,5,6,7,8,9,10}. 
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Actually, 

8({1, 2, 5, 6, 10}) = {1,2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14, 15, 20, 21, 24, 27, 32, 35, 

44, 48, 49, 54, 55, 63, 80, 98, 99, 120, 125, 175, 224, 
242, 384, 440, 539, 2400, 3024, 4374, 9800} 
contains all of the solutions to ([[]) with y = 7. It is easily checked that 

2 /r,\ -2 /r \ -4 



4374 
4375 



so that the fraction yy, which corresponds to b = 10, does not appear; but in finding the 
above decomposition by tracking the iterative procedure in computing 6({1, 2, 5, 6, 10}) this 
fraction made an appearance and then is later canceled out. 

Looking at Lehmer's tables we find for all primes p up to 41, with the exception of p = 7 
and p — 41 that 

z p :=S({l,2,...,p}) 

coincides exactly with the complete set of p-smooth solutions of ([1]). In the case p = 41, we 
found that Z41 has 890 of Lehmer's solutions; it is missing the largest solution, 63927525375. 
Note that 

63927525375 3 3 x 5 3 x 7 7 x 23 



The least n for which 



63927525376 2 13 x ll 4 x 13 x 41' 
63927525375 E 5{{1, 2, . . . , n}) 



is n = 52. 

We calculated Zigg in a week on a PC using Mathematica. It has 346192 elements. A 
histogram of the logarithms of these numbers is shown in Figure 1. 

The prime before 199 is 197. The number of new b in z 19 g, i.e. those not in z 197 , is 43215. 
For all but 300 of these b it is the case that 197 | b{b + 1). The other 300 have a smaller 
largest prime factor. Here is a table of the number of these b sorted by the largest prime 
factor of b{b + 1): 



Prime 


127 


131 


137 


139 


149 


151 


157 


163 


167 


173 


179 


181 


191 


193 


Number 


1 


1 


2 


2 


5 


5 


10 


8 


13 


16 


25 


36 


43 


51 



197 
82 



3. Solving the Diophantine equation 

In this section we analyze the solutions of 

b B + l /3 
6+1 X B ~ /3 + 1 

This equation can be written as 

6(fl + l)(/9 + l) = {b+l)B(3. 
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Figure 1. A histogram of the set {log b : p | b(b + 1) 
with the PDF of the normal distribution. 
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Figure 2. A plot of (n, \z Pn \ — |-2 P „_J) where p n denotes the nth prime. 



Let gcd(6, B) = g and put 

b = gu B = gv 

where clearly gcd(w,f) = 1. Then we have 

u(B + l)(/3 + l) = (b + l)v(3. 
Now, gcd(w, (6 + l)v) = 1 so it must be the case that u \ {3 , say (3 = hu. Then we have 

(5 + 1)03 + 1) = (b + l)vh. 
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Then v \ (f3 + 1), say (3 + 1 = vx and h \ (B + 1), say B + 1 = /«/. Thus, in terms of the 
variables g, h, u, v, x, y, we have the three equations 

hy - gv =1 
vx — hu =1 
xy — gu =1 

along with b = gu, B = gv and (3 = hu. If we solve for x and y in the first two equations 
and substitute the results into the third equation we have 

uh + gv = hv — 1. (3) 

Notice that if h and v are given positive coprime integers then there exist unique positive 
integers u and g that satisfy ([3]). Then x and y can be determined as well as b,B and (3. 
Thus, a coprime pair (h, v) leads to a triple (6, B, (3) and vice- versa. 

From another perspective, one may ask how to find a pair (6, B) from a given (3. Take any 
divisor w of (3 and any divisor v of /3 + 1. Then /i = and x = (/3 + l)/f, from which we 
see that 

v — u h — x 



y = T~ 9 



vx — hu vx — hu 

Thus, the pair (u, v) leads to the pair (b, B) given by 

uh — ux n u . n , . „ hv — xv v n . n „ . f 6 

6= — =(3 --((3 + 1) B = — = -13- {(3 + 1 ) = -, 

vx — hu v vx — hu u u 

where we assume that u < v. And we've already seen that given b < B we can define 

b B 



u 



(b,B) (b, B) ' 

Thus, given (3, there is a one-to-one correspondence between pairs (u, v) with u < v, u \ (3, 
v | (/3 + 1) and pairs (6, 5) with b < B and x = -^j. 

4. Trees 

In this section we illustrate the way in which 601425 enters into Z20 = <5({1, 2, . . . 20}). 
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This tree illustrates the factorization 
601425 



601426 




6\ [7\ 1 /8\ Via /i5\ fn\ 1 /is\ /19 



7J\8J \9 J \12J\1QJ\18J V 19 7 V 20 

5. LUCA AND NAJMAN'S WORK 



Luca and Najman published a paper in Mathematics of Computation in 2011 which im- 
proved on Lehmer's previous work [LN]. They gave a list of 13325 pairs of numbers which 
were 97-smooth neighbors and claimed that this list was complete. Their computation took 
15 days on a computer with a 2.5 GHz processor. By comparison, when we searched for 
97-smooth neighbors, our method took just under 20 minutes and, much to our surprise, 
produced 13333 neighbors. In comparing their list with ours, Luca and Najman found 37 
solutions not on our list. However, we found 45 solutions not on their list! Independently of 
our work, Luca and Najman found the 49 solutions they had missed previously. 

Our calculations of 199-smooth numbers took more than a week and produced 345192 
solutions, some of which were the solutions we had missed in our 97-smooth calculation. In 
fact, our ^199 is only missing one solution from the Luca-Najman list, 9591468737351909375. 

6. The largest Luca-Najman solution 

The 79-smooth neighbor pair which starts with 

(3 = 9591468737351909375 

was found by Luca and Najman but is not on our list (which we call z 199 ) of 199-smooth 
neighbors. We can prove that it will first appear in our method when we search for 227- 
smooth neighbors. 

Given a number (3, there is a 1-1 correspondence between pairs (b, B) with b(B + l)/((6 + 
l)B) = (3/((3 + 1) and pairs (u,v) with u \ (3 and v \ ((3 + 1). In one direction this 
correspondence is given by 

a 

b = 0--O3 + 1) 
B = ^-0 + 1) = !* 

U U 

There are 1440 divisors of (3 and 5632 divisors of (3 + 1. This means that there are a total 
of 8110080 pairs (u,v) to consider. For each of these pairs we computed the pair (b,B) and 
then computed the largest prime factor of b(b + 1)B(B + 1). The least out of all of these 
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largest prime factors was p = 227 which appeared for several pairs (b,B), in particular for 
the pair 

b = 285406166331883519 B = 294159243066390624. 

Therefore, /3 cannot appear as a solution in our method unless we go up to p = 227. 

Further analysis led us to find an explicit tree for /3 all of whose bottom nodes are either 
in ^199 or else appear within the first two iterations (which are quick to compute) arising in 
the computation of z 2 27', thus we can show that f3 G ^227 without computing all of ^227. 

We now describe the tree in detail. Given a (3 we have described a 1-1 correspondence 
between pairs (u,v) and pairs (b,B). For a pair (b,B) let us use the notation 

F(b,B)=P 

to denote that 

b (£ + !) _ g 
(6+1) X B ~(/3 + l)' 

Thus, 

F(285406166331883519, 294159243066390624) = 9591468737351909375. 
is the first step of our computation. Let's use 

b = 285406166331883519; B = 294159243066390624; (3 = 9591468737351909375 
for the rest of this section, so that F(b, B) = (5. Now let 

gi = 2229716045541599; g 2 = 2247272709023744 



h x = 186642247267999999; h 2 = 510640590102749183. 
F( gi ,g 2 ) = b; F(h 1 ,h 2 ) = B. 



and 
Then 

Moreover, 

9l £ ^199 h 2 G ^199 

so that we may restrict our attention now to g 2 and hi. Next, we let 
i t = 907177810312319; i 2 = 911608699868750 

and 

ji = 1671690051584; j 2 = 1672934505788. 

Then 

F(i 1 ,i 2 ) = h 1 F(j 1 ,j 2 )=g 2 
and j 2 G ,2199; we are left to decompose ii,i 2 , and j\. Let 

h = 341611712; k 2 = 341681535; m 1 = 300775; m 2 = 301040. 
These satisfy 

F(ki, k 2 ) = ji F(mi, m 2 ) = k 2 
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and ki,m 2 G £199. Let 

ni = 3405; n 2 = 3444; o x = 454; o 2 = 524; 

then 

F(n 1 ,n 2 ) = mi, F(o 1 ,o 2 )=n 1 

and n 2 ,o 2 G ,2199. Let 

Pi = 199802399641; p 2 = 199846415040; q x = 16503580; q 2 = 16504943; 

then F(p 1 ,p 2 ) = i ± ] F(q 1 ,q 2 ) = p 2 ; and p x G 2:199. Let 

n = 561824; r 2 = 581624; Sl = 49664; s 2 = 54480; 

then F(r 1 ,r 2 ) = q ± ; F(si,s 2 ) = r±; and r 2 , Si G 2199. Let 

h = 6810; t 2 = 7783; u x = 108033083250000; u 2 = 122557101693480; 

then F(t 1 ,t 2 ) = s 2 ; F(u l7 u 2 ) = i 2 ; and t 2 ,Ui G z 199 . Let 

v l = 10638314820; v 2 = 10639238337; w 1 = 1451240; w 2 = 1451438 

then F(vi,v 2 ) = u 2 \ F(wi,w 2 ) = v\\ and v 2l w\ G Z199. Let 

x 1 = 92852; x 2 = 99198; y x = 4312; y 2 = 4508; 

then F(x 1 ,x 2 ) = w 2 ; F(y l ,y 2 ) = x 2 ; and x X) y 2 G z 199 . Finally, o 1 ,t 1 ,y 1 G z 227 . This last 
fact is determined by computing the first two iterations in the process for determining 2227 
which only takes a few seconds. 

7. Computational time 

We have computed 346192 solutions to 

p I 6(6 + 1) p < 199. 

199 is the 46th prime. Would it be possible to use the Pell's equation method to find the 
complete list of 199-smooth neighbors? No. To use the Pell's equation method [LI] to find 
all of the solutions, one would first have to find the fundamental solutions of 

2 46 - 1 = 70368744177663 

different Pell's equations and then check up to the 100th solution of each. The difficulty 
in finding the continued fraction expansion of \fd grows dramatically with d. We expect 
that for the prime 199 we would find d such that the period of \fd is as large as 10 40 . Even 
using the subexponential algorithm employed by Luca and Najman [LN] it would still be 
impossible to go this far. 

Would it be possible by this method to find all of the solutions we found? No. For the 
numbers b on our list £199 the continued fraction expansions of \/b(b + 1) have unusually 
small periods; in fact the largest such period is only 38. Suppose one does a calculation where 
they check all 2 46 square-free gPs that are 200-smooth to see which lead to short continued 
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fraction expansions as described above; even this calculation we estimate would take, at a 
conservative minimum, at least 3000 years on the computer we used. 

8. Largest solutions 

Here is a list of (q, 6) where q is a prime number up to 197 and 6 is the largest element of 
Z199 for which q | (6(6 + 1) and p | 6(6 + 1) =>- p < q. 

(2, 1), (3, 8), (5, 80), (7, 4374), (11, 9800), (13, 123200), (17, 336140), (19, 11859210), 

(23, 5142500), (29, 177182720), (31, 1611308699), (37, 3463199999), (41, 63927525375), 

(43, 421138799639), (47, 1109496723125), (53, 1453579866024), (59, 20628591204480), 

(61, 31887350832896), (67, 12820120234375), (71, 119089041053696), 

(73, 2286831727304144), (79, 1383713998733898), (83, 17451620110781856), 

(89, 166055401586083680), (97, 49956990469100000), (101,4108258965739505499), 

(103, 19316158377073923834000), (107, 386539843111191224), 

(109, 90550606380841216610), (113, 205142063213188103639), 

(127, 20978372743774437375), (131, 1043073004436787852800), 

(137, 65244360004072055000), (139, 152295745769656587384), 

(149, 6025407960052311234299), (151, 1801131756071318295624), 

(157, 277765695034772262487), (163, 1149394259345749379424), 

(167, 2201197005772848768608), (173, 4574658033790609920000), 

(179, 9021820053747825025975), (181, 13989960217958128903124), 

(191, 75121996591287627735039), (193, 444171063468653314858175), 

(197, 25450316056074220028640), (199, 589864439608716991201560) 

Note that the Luca-Najman number 9591468737351909375 would correspond to q = 79 
above, where we have 1383713998733898. 

9. Smoothness and the ABC-conjecture 
The ABC equation is 

A + B = C. (4) 

The 1985 ABC-conjecture of Masser and Oesterle asserts that for any e > there is a 
k(c) > such that for all solutions to (j4]) the inequality 

C < /?(e)rad(ABC) 1+e 

holds, where rad(n) := n p |nP> caue d the radical of n, is the product of the prime divisors 
of n. In studies of the ABC conjecture solutions to (HI) there are a variety of indicators used 
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to measure the solutions. A chief one is the quality: 

logC 



q(A,B,C) 



lograd(ABC)' 

In fact, it seems to be standard to call any triple (A, B, C) whose quality is at least 1 "an 
ABC-triple." Much work has been done to systematically find all ABC-triples. The first 
22763667 ABC-triples are available for download from the abc@home web-site [A]. These 
include all ABC-triples with C < 10 18 . 

Lagarias and Soundararajan [LS] define what they call the smoothness exponent k of an 
ABC solution by 

k {A,B,C) = - — - — - 
log log C 

where z is the largest prime factor of ABC, and they define the number kq by 

k,q := liminf Kq(A, B, C). 

C— >oo 

They observe that the ABC conjecture implies that k,q > 1 and conjecture that k = 3/2. 
They also prove that the Generalized Riemann Hypothesis implies that Kq < 8. 

We introduce another measure of the size of a z-smooth solution of with a quantity 
we call the smoothness index given by 

S( A,B,C).J^ 
log z 

where we take z to be the largest prime factor of ABC . For example, 

2 + 25 = 27 

has C = 27 and the largest prime factor of 2 x 25 x 27 = 2 x 3 3 x 5 2 is 5 so that s(2, 25, 27) = 
M = 2.04782. 

log 5 

Remark 1. For comparison, it is easy to make the smoothness exponent large (just take C 
to be a large prime ) but hard to make it small. For the smoothness index, the challenge is to 
make it large. 

Remark 2. The quality is more a measure of the average size of the prime divisors of ABC 
whereas the smoothness index measures the size of the maximum prime divisors of ABC . 
Note that an ABC-triple can have a high quality even if a few of the prime divisors are large 
whereas that cannot happen for measures of smoothness. 

Remark 3. The ABC conjecture asserts that q — > 1 as C — > oo. However, by work of Balog 
and Wooley [BW] s(A,B,C) is unbounded. 

Remark 4. Work of Lagarias and Soundararajan [LS], leads us to expect triples with 
s(A,B,C) about as large 

Thus, theoretically we know s(A,B,C) can be arbitrarily large; for us the challenge is 
actually finding (A, B, C) with s(A, B, C) large. In other words it's a computational challenge 
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rather than a theoretical challenge to find large values of s(A,B,C). Our algorithm finds 
large values of 

s(l,B,B + l). 

The largest value we found is for 

B = 19316158377073923834000 

we have 

s(l,B,B + l) = 11.0719 

The value of B here is ~ 1.9 x 10 22 which is beyond the range where systematic study of the 
ABC conjecture has been conducted. We know of no larger values of s. 

10. Some data about smoothness 

We call a triple (A,B,C) maximally smooth if s{A 1 ,Bi,C\) < s(A,B,C) for all C± < C 
(or for Ci — C and A\ < A). Here are the first thirteen maximally smooth triples, together 
with their smoothness index: 



1 + 3 


= 4 


1.262 


3 + 5 


= 8 


1.292 


1 + 8 


= 9 


2 


2 + 25 


= 27 


2.048 


5 + 27 


= 32 


2.153 


1 + 80 


= 81 


2.730 


3 + 125 


= 128 


3.015 


32 + 343 


= 375 


3.046 


49 + 576 


= 625 


3.308 


5 + 1024 


= 1029 


3.565 


1 + 2400 


= 2401 


4 


1 + 4374 


= 4375 


4.308 


7168 + 78125 


= 85293 


4.427 



The triples on this list seem very appealing. In particular, 

3 + 5 3 = 2 7 

1 + 5 x2 4 = 3 4 

2 5 + 7 3 = 3 x 5 3 

5 + 2 10 = 3 x 7 3 

1 + 3 x 2 7 = 7 4 

7 x 2 10 + 5 7 = 13 x 3 8 
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are all of the form 

ax 1 + by m = cz n 

with a, b, c, x, y, z each being 1 or a prime and £,m,n > 2. 

A triple (A, B, C) is called an "ABC-triple" provided that C > rad(ABC). The first 
22763667 triples are available for download from the abc@home web-site [A]. These include 
all of the triples for C < 10 18 . The largest C for a triple on this list is for the triple 

131854322526743011 + 9091517323167918864 = 9223371645694661875 

or in factored form 

13 4 x 16651 3 + 2 4 x 3 8 x 53 4 x 3313 2 = 5 4 x 7 7 x 37 x 59 2 x 373 2 . 

This triple has smoothness index 4.493. Here we have C = 9.2 x 10 18 . 

Among these more than 22 million ABC - triples the triple with largest smoothness index 

is 

176202799695875 + 3178472661789594624 = 3178648864589290499 
with smoothness index 11.0653. In factored form this is 

(5 3 x 23 3 x 42 5 ) + (2 10 x 3 7 x 7 6 x 13 3 x 17 2 x 19) = ll 4 x 31 x 37 4 x 43 3 x 47. 



11. Lehmer's table 
For each prime z < 41 and each k < 6 Lehmer [L2] found the largest n such that 

p | n(n + 1) . . . (n + h — 1) =>- p < z. 

For example, for z = 37 and h = 3 he found n = 17575 which means that 17575 x 17576 x 
17577 = 2 3 x 3 4 x 5 2 x 7 x 13 3 x 19 x 31 x 37 has no prime factor larger than 37 and that 
17575 is the largest number with that property. Our algorithm produces many more such 
data points, but in each case we can only say that our number is a lower bound for the actual 
number, since we don't know for any z > 41 whether our list of b is complete. 
But to give some examples, we found that 

£(134848 x 134849 x 134850) = 43 
£(192510 x 192512 x 192512) = 47 



£(1205644 x 1205645 x 1205646) = 53 
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Below is a table of lower bounds. 
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oi oqcn 
ZlZooU 


lOoU 


i ci 7 
101 / 


Ofi K. 

Zoo 




/I Q 

4o 


/lOI 1 QS7QQKQQ 

4ziioo / yyooy 


oi oqcn 
ZlZooU 


lUo /o 


i ci 7 
101 / 


Q/i n 
o4U 




A 7 
4 / 


1 1 no/iQK79Qi oc; 
iiuy4yo / zoizo 


ZlZooU 


1/0/0 


i ci 7 
101 / 


Q/i n 
o4U 


1 QA 
lo4 


53 


1453579866024 


1205644 


17575 


1517 


340 


184 


59 


20628591204480 


1205644 


17575 


1517 


528 


527 


61 


31887350832896 


1205644 


17575 


1767 


528 


527 


67 


31887350832896 


2018978 


17575 


5828 


1271 


527 


71 


119089041053696 


3939648 


17575 


5828 


1271 


527 


73 


2286831727304144 


3939648 


70224 


5828 


1271 


527 


79 


2286831727304144 


15473808 


70224 


5828 


3476 


527 


83 


17451620110781856 


15473808 


97524 


5828 


4897 


4896 


89 


166055401586083680 


407498958 


97524 


5828 


4897 


4896 


97 


166055401586083680 


407498958 


97524 


7565 


7564 


4896 


199 


589864439608716991201560 


768026327418 


61011223 


1448540 


44250 


18904 



We note some minor errors in Lehmer's tables: his h = 5 entry for p = 41 should be 1517. 
Also, his Table IIA is missing the entries 10935, 12901781, and 26578125 in the 29 column 
and 4807, 12495, 16337, 89375 from the 31 column. 

12. Future directions 

The unanswered question in all of this is why does it work? Can we really get all of the 
b for which b(b + 1) has only small prime factors by repeatedly applying 5 to a small initial 

set? 

A reasonable thing to do would be to run this for a year on a supercomputer and generate 
millions of smooth neighbors, perhaps some with as many as 30 digits. 

It is true that computations get longer as you raise your smoothness barrier. But the length 
of computation time is due mainly to the vast number of solutions, not the complexity in 
finding them. 
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The same technique we have elucidated above works surprisingly well to find solutions to 
the more general problem 

p \ b(b + k) =>- p < z 
when k is an odd integer. Fix a difference k. Start with an initial set S and let 

then S" = (S')' and so on. Iterate this procedure until it stabilizes and call the result 

5 k (S). 

It would be interesting to do extensive computations with a range of differences k. 

For even integers the process has to be modified somewhat. For example, to generate 
solutions of 

p | b{b + 2) =>- p < z 

one should start with the set z p = <5i({l, 2, . . . ,p}) for some p and then look for solutions of 

b B + l g 
6+1 X B ~ /3 + 2 

with b,B G z p . 
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